;***************************************************************
;discription： ndvi最小二乘法方程中的slope斜率分析程序入点
;采用最小二乘法多项式拟合，函数公式为：y=a+bx;slope就是b参数
;本程序完成了最小二乘法方程拟合运算工作
;
;author: wangqiang cdut 2016-10-1
;email：  cqwensten@163.com  419436974@qq.com
;***************************************************************
pro pro_leastsquare_slope
  ;严格编译规则
  COMPILE_OPT IDL2
  ;二次开发模式初始化，恢复envi的sav文件
  ENVI,/RESTORE_BASE_SAVE_FILES
  ;初始化envi，保存日志文件
  ENVI_BATCH_INIT,log_file='ndvi_slope.log'
  starttime=systime(1)
  PRINT , '-----' ,starttime, '-----'
  ;********************************开始****************************************

  ndvi_file= envi_pickfile(title='选择NDVI时间序列文件');这个函数可进行文件选择
  ;open the image file
  envi_open_file, ndvi_file, r_fid=ndvi_fid ; open the ndvi file
  if (ndvi_fid eq -1) then begin
    envi_batch_exit
    return
  endif
  ;获取行、列号，波段数，数据类型
  envi_file_query, ndvi_fid, ns=ncols, nl=nrows, nb=nbands, Data_Type=data_type, xstart=xstart, ystart=ystart
  print,'the ndvi file info（FID,列数，行数，波段数） :',ndvi_fid, ncols, nrows, nbands

  ;获取文件投影
  map_info = envi_get_map_info(fid=ndvi_fid)
  help,map_info
  print,map_info
  print,ndvi_file

  ndvi_data = FltArr(ncols,nrows,nbands)    ;  FLOAT  Floating point
  dims = [-1, 0,ncols - 1 , 0,nrows - 1]
  FOR i=0,nbands-1 DO BEGIN
    ;get the image data in the openning file which band pos=i
    ndvi= envi_get_data(fid=ndvi_fid, dims=dims, pos=i) ;波段数从0开始,提取出训练数据,并且第一层数据为目标数据
    ndvi_data[*,*,i] = float(ndvi)
  ENDFOR


  fitDegree = 1
  xlist = indgen(nbands) + 1;
  slope_data = make_array(ncols,nrows,type=5,value=0.0)

  ;采用最小二乘法多项式拟合，函数公式为：y=a+bx;
  FOR yindex=0, nrows-1 DO BEGIN
    FOR xindex=0, ncols-1 DO BEGIN
      ylist = ndvi_data[xindex,yindex, *]
      slope = 0.0
      
      idx_complement = WHERE(ylist gt -30.0, count, complement=idx_ncomplement, ncomplement=ncount)
      ;无数据区
      if ncount gt 0 then begin
        ;slope = -9999.0
        slope = -1
      endif else begin
        slope = FuncLineSlope_Cal(xlist,ylist)
      endelse

      ;保存结果
      slope_data[xindex,yindex] = slope
    ENDFOR
  ENDFOR

  ;保存模型结果
  output_file = DIALOG_PICKFILE(TITLE = 'ndvi slope输出结果...')
  ;ENVI_WRITE_ENVI_FILE, npp, OUT_NAME = output_file,Map_Info = map_info
  ;二进制方式输出
  IF strlen(strtrim(output_file,2)) GT 0 THEN BEGIN
    OPENW,lun, output_file,/get_lun
    ;写入npp结果数据到文件中
    WRITEU,lun,slope_data
    FREE_LUN,lun
    ;写出头文件信息
    ENVI_SETUP_HEAD,fname=output_file, ns=ncols, nl=nrows, nb=1, interleave=0, data_type=5, offset=0, map_info=map_info, /write, /open
  ENDIF
  ;********************************退出***********************************
  ;print the end time
  endtime=systime(1)
  duration=endtime-starttime
  PRINT ,'-----' ,endtime, '-----'
  ;关闭envi二次开发模式
  ENVI_BATCH_EXIT
end

;*********************************************
;
;计算像元的斜率
;*********************************************
FUNCTION FuncLineSlope_Cal, X, Y
  slope = 0.0
  ;最小二乘法计算直线斜率
  results = poly_fit(X,Y,1,/DOUBLE,STATUS=status)
  IF(status eq 0) THEN BEGIN
    ;y=a+bx,其中b就是斜率，a为截距
    ;a=results[0];b==results[1]
    slope = results[1]
  ENDIF

  ;返回斜率
  return,slope
END
